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We theoretically investigate Coulomb drag in a system of two parallel monolayers of graphene. 
Using a Boltzmann equation approach we study a variety of limits ranging from the non-degenerate 
interaction dominated limit close to charge neutrality all the way to the Fermi liquid regime. In the 
non-degenerate limit we find that the presence of the passive layer can largely influence the conduc- 
tivity of the active layer despite the absence of drag. This induces a non-trivial temperature behavior 
of the single layer conductivity and furthermore suggests a promising strategy towards increasing 
the role of inelastic scattering in future experiments. For small but finite chemical potential we find 
that the drag resistivity varies substantially as a function of the ratio of inelastic and elastic scatter- 
ing. Furthermore, we explicitly show that the clean system has a well-defined drag resistivity even 
though the individual conductivities diverge. We find that an extrapolation from finite chemical po- 
tential to zero chemical potential and to the clean system is delicate and the order of limits matters. 
While the drag resistivity pd extrapolates to zero upon taking the limit lima-,.^ lim Ma=|[lp _ > .o Pd = 
it has a finite value in the opposite order of limits lim Mo=Pj) _ > .o lrnia-»oo pd = — ^ (fJ-a and p p are 
chemical potentials of the active and passive layer). The limiting value in the latter case is set 
by the interaction dominated single layer conductivity ao of clean graphene and in that sense is a 
universal number. In the Fermi liquid regime we analyze drag as a function of temperature T and 
the distance d between the layers and compare our results to existing theoretical and experimental 
results. In addition to the conventional l/d 4 -dependence with an associated T 2 -behavior we find 
there is another regime of l/d 5 -dependence where drag varies in linear-in-T fashion. The relevant 
parameter separating these two regimes is given by d — Td/vF (vf is the Fermi velocity), where 

<gC 1 corresponds to T 2 -behavior, while d> 1 corresponds to T-behavior. We speculate that the 
broad crossover between these two regimes was observed in recent experiments on graphene as well 
as old experiments on conventional two dimensional electron gases. We close with a discussion of 
the role of screening and the determination of the drag resistivity as function of the charge carrier 
densities in the two layers under very general circumstances covering the whole crossover from the 
non-degenerate to the degenerate limit in both layers independently. 



I. INTRODUCTION 

Graphene, a two dimensional system of carbon atoms 
arranged on a hexagonal lattice with an emerging Dirac 
type low-energy electronic dispersion continues to attract 
considerable interest on the theoretical and experimental 
front 1 . One remarkable feature in experiments is that 
so far they have revealed only very limited information 
about interactions. The most prominent manifestations 
of interaction effects arc the observation of the fractional 
quantum Hall effect 2 ' 3 as well as the logarithmic scaling 
of the Fermi velocity of the Dirac particles which was 
recently seen in quantum oscillation measurements on 
ultra-clean suspended samples 4 ' 5 . However, with ever in- 
creasing sample quality one expects to eventually be able 
to reach the hydrodynamic collision-dominated regime 6-9 
allowing to observe non-trivial many-body physics such 
as a collective cyclotron resonance 10 ' 11 or an anomalously 
low viscosity 12 . Also, a quantum-critical version of the 
Kondo effect possibly comes within reach 13 . A very di- 
rect manifestation of Coulomb interactions is provided by 
Coulomb drag experiments, the effect of electrons moving 
in one plane dragging along electrons in a plane parallel 
to the one in which the current is driven. This effect 
has a long history in the context of two dimensional elec- 
tron gases 14-24 . In graphene this problem has previously 



been studied in experiment 26-28 and in a recent series of 
theoretical works 29-37 . Here we report on theoretical re- 
sults in the framework of a Boltzmann approach. Our 
approach goes beyond former approaches in that we al- 
low for varying single layer properties as a function of all 
parameters. While this is not vital in the description of 
Coulomb drag in the Fermi liquid regime, \n/T\ ^> 1, this 
becomes crucial in the non-degenerate limit, |/x/T| <C 1, 
where interaction effects can dominate the single layer 
properties and an interesting interplay between elastic 
and inelastic scattering can be observed. Experimen- 
tally, there are indications that this regime should be 
within reach in experiments using samples prepared on 
hexagonal boron nitrid substrates where due to the atom- 
ically smooth surface that is relatively free of dangling 
bonds and charge traps high purity can be achieved and 
the puddle regime can be suppressed to very low densi- 
ties 38 . Our approach, like other theoretical approaches to 
date, is not valid in this strongly inhomogeneous regime. 
Throughout the paper we keep our results as general as 
possible, meaning we try to keep the number of Dirac 
cones N in final expressions, if possible. This implies 
that our results should equally apply to three dimensional 
topological insulators, whose surfaces are characterized 
by an odd number of Dirac cones (in the case of weak 
topological insulators there is an even number of Dirac 
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cones). A possible drag setup in such a system is even 
more straightforward and very natural in the sense that 
slab systems with a finite size gap in z-direction host a 
natural setting in which our results apply. However, we 
stress that the localization physics in these theories is 
different due to the helical nature of the surface Dirac 
fermions. 



A. General properties of drag 

In the experimental setup, Fig. 1, two monolayers of 
graphene are separated by a distance d. We assume 
that in-between the monolayers there is an insulating 
region filled with a dielectric with a dielectric constant 
e r . Throughout the paper the dielectric constant e r is 
not a function of the vertical coordinate. This situation 
has been studied elsewhere 35 . We assume that the two 
layers can be individually gated such that the carrier con- 
centration in both layers can be adjusted independently. 
Furthermore, we divide the two layers into active and 
passive layer, where active layer refers to the fact that 
within this layer a current is driven, while the passive 
layer will not carry current. In a standard experiment a 
current I\ is driven through the active layer. If no cur- 




FIG. 1: Schematic setup of a drag experiment. In the active 
layer (a) a current Ii is driven. In the passive layer (p) a 
voltage V2 is induced such that overall there is no current 
flow in the passive layer. The drag resistance is defined as 
Ra = -V a /h. 

rent is allowed to flow in the passive layer this induces a 
voltage difference V2, allowing to define a drag resistance 
R 2 = -V 2 /h. 

We consider the response tensor which has a structure 
similar to the one in the Hall effect. We apply the electric 
field E a only in the active layer a but induce a current 
in the active layer a called j a as well as j p in the passive 
layer p. Consequently, there are layer-diagonal and off- 
diagonal conductivities involved: 

(t)- (2 ;')•(•)■ <" 



which includes the individual conductivities a a of the ac- 
tive and that of the passive layer, er p , while the drag 
conductivity is denoted a 4. In the concrete experiment, 
however, the boundary conditions are different and the 
passive layer does not carry current. Demanding j p = 
requires a field counteracting the flow in the passive layer 
which is given by E p = — ^-E a . This implies that the 

drag resistivity (or transresisitivity) is given by 



Ja 'top cr a cj p — a d 

It is important to realize that like in the case of thermal 
transport pd can be finite even if the individual conduc- 
tivities er a , <7 P , and aa diverge, which we show explicitly. 
This is an effect of the boundary condition of vanishing 
charge flow in the passive layer analogous to a finite ther- 
mal conductivity in thermal transport in Fermi liquids. 



B. Summary of results 

Graphene bilayers turn out to provide an exceptionally 
versatile arena in which one can theoretically as well as 
experimentally vary a large number of parameters: (i) 
the temperature T, (ii) the chemical potentials of the in- 
dividual layers fj, a and fx p , (hi) the interlayer spacing d, as 
well as (iv) disorder and (v) interaction strength via the 
dielectric environment. Within this work we do not at- 
tempt to exhaust all the possibilities offered by the above 
parameters but instead concentrate on the most interest- 
ing situations. Our main results concern among others 
the interplay of interactions and disorder in the limit of 
small chemical potentials, the so-called non-degenerate 
limit where |/z a | , \fi p \ <C T in both layers individually. 
We furthermore describe the crossover from this inter- 
esting limit to the more conventional Fermi liquid regime 

On a different note we study the dependence on dis- 
tance d in some detail in both the non-degenerate and 
the degenerate limit. In experiments distances d can be 
realized for which the regime d -C vf/T can be reached 
in absence of leakage currents for reasonable tempera- 
tures. This is an interesting limit especially in the non- 
degenerate case since the typical momentum of electrons 
involved in scattering events is on the order T/vf and 
consequently the interlayer Coulomb interaction can be 
considered essentially undamped, see Eq. (18), meaning 
we have the maximal effect of the inter-layer interac- 
tion strength and inelastic scattering can effectively be 
increased. However, we also study the opposite limit, 
d 3> vp/T, which is particularly interesting in the Fermi 
liquid regime revealing a formerly not discussed regime. 

Our main findings are as follows: (i) For zero doping in 
both layers, i.e., chemical potential fi a .p = 0, the passive 
layer remains in equilibrium and consequently the drag 
resistance is zero. This happens by virtue of the current 
carrying state being of zero total momentum which is 
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enforced by particle-hole symmetry. Despite the passive 
layer remaining in equilibrium we find that its presence 
still has a large effect on the active layer. In the limit 
Td/vF <C 1, the conductivity of the active layer is re- 
duced by roughly a factor of two. This comes about 
due to scattering of the electrons and holes in the active 
layer from plasmons in the passive layer which becomes 
increasingly pronounced upon decreasing temperature. 
The effect is unusually large since in a normal Fermi liq- 
uid the influence of the passive layer on the transport 
properties of the active layer is negligibly small. Overall, 
we find that the active layer shows a non-trivial tem- 
perature dependence of the conductivity. This provides 
a potential route towards increasing inelastic scattering 
in the active layer in a setup where the active layer is 
sandwiched between a (possibly large) number of pas- 
sive layers which act as a reservoir for inelastic scatter- 
ing, (ii) In the non-degenerate limit \[i a \, |MpI ^ T 
we find an interesting crossover from disorder dominated 
drag to interaction dominated drag. In the limit of zero 
disorder, in which all individual conductivities diverge, 
the drag resistivity remains finite, which we show explic- 
itly. Furthermore, we find that the approach to the clean 
limit combined with the approach of the Dirac point is 
subtle and the order of limits matters: this means that 
taking disorder to zero and subsequently the chemical 
potential to zero or the other way around yields dif- 
fering results. In the former case one obtains a finite 
drag resistivity even at the Dirac point which is given 
by the inverse single layer conductivity of the clean sys- 
tem. In the latter order of limits one obtains zero, (iii) 
In the Fermi liquid regime we study the dependence of 
the transresistivity on the distance d between the two 
layers. We find that there is an interesting crossover in 
the behavior of the drag resistivity as a function of the 
distance d which has previously not been discussed in 
the context of Fermi liquids 18,19 , namely a change from 
l/d 4 -behavior to 1/d 5 -behavior. This effect is accom- 
panied by a crossover in the temperature dependence 
which goes from T 2 to T, which could be relevant for the 
understanding of recent experiments in the Fermi liquid 
regime as well as older experiments on conventional two- 
dimensional electron gases. As a byproduct we derive the 
standard formula for Coulomb drag in the Fermi liquid 
regime 18 ' 19 from a very simple one-mode approach to the 
Boltzmann equation, which to the best of our knowledge 
has not been discussed before, (iv) We describe the full 
crossover from the non-degenerate limit into the Fermi 
liquid regime |/x a |>|Mpl ^ T in which the single layer 
conductivities are disorder dominated. We find that the 
effect of screening in this limit brings us towards orders 
of magnitude of the drag resistivity which are very com- 
patible with experimental results. 

Technically we use the kinetic approach, which re- 
quires the full numerical solution of coupled Boltzmann 
equations for the distribution functions of the electrons 
and holes in the individual layers, thus for four coupled 
Boltzmann equations. This is a straightforward but non- 



standard application of the variational principle 39 and 
consequently explained in some detail. The present work 
goes beyond former theoretical works in mainly three as- 
pects: (a) We do not use a relaxation time approximation 
but instead solve the Boltzmann equation numerically 
within a two-mode approach (this is logarithmically ex- 
act in the strong coupling limit), (b) The description 
of the interplay of interactions and disorder especially in 
the non-degenerate limit is facilitated by the two-mode 
description, which is the minimal number of modes re- 
quired for a faithful account, (c) We do not take the 
individual layer conductivities as input parameters but 
instead calculate them for every set of parameters which 
leads to qualitative and quantitative changes in the non- 
degenerate limit. 



C. Organization of the paper 

The organization of the paper is such that we start 
with a discussion of the setup in Sec. II , which includes 
a discussion of the model Hamiltonian (Sec. II A). In 
Sec. II B we first discuss the sources of current relax- 
ation in Sec. II Bl, then the associated time scales in 
Sec. II B 2, as well as the effect of screening in Sec. II B 3. 
The generic framework of the Boltzmann equation and 
the matrix formalism used follows in Sec. III. We first 
introduce the coupled kinetic equations necessary to de- 
scribe drag in Sec. Ill A. Here we also explain how the 
effect of drag manifests itself in the structure of the cou- 
pled equations in linear response. We then move towards 
the variational ansatz in Sec. MB and shortly review the 
variational principle. We also discuss the minimal num- 
ber of modes required for a faithful description within our 
problem. In a last step, Sec. Ill C, we introduce a generic 
matrix formalism derived from the variational principle 
which enables us to calculate drag from an inversion of a 
matrix. Readers not interested in technical details may 
skip Sec. Ill and directly move to the results: we start 
with a discussion of the non-degenerate limit in Sec. IV. 
In a first step we discuss the case of both layers at the 
Dirac point in Sec. IV A. The results lead us to propose 
an experimental setup, in which the effect of inelastic 
scattering can be increased considerably in Sec. IV B. We 
then move to finite chemical potential, but still T » fi, 
in Sec. IV C. There we discuss that the transresistiv- 
ity can be finite even in a clean limit and show that drag 
largely depends upon the ratio of elastic to inelastic scat- 
tering. We furthermore discuss that the extrapolation to 
zero density and zero impurity density is delicate and 
the order of limits matters. In a next step we analyze 
the Fermi liquid regime in Sec. V which is what has been 
analyzed predominantly in other works. We first give an 
alternative derivation of the standard results of drag in 
Fermi liquid theory 18 ' 19,34 ' 35 obtained in the Kubo ap- 
proach in terms of coupled Boltzmann equations. Then 
we discuss a number of different situations and study the 
behavior of drag with distance d in great detail. Im- 
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portantly, we find a regime of temperature and distance 
dependence which has been overlooked in previous works 
and could be important for the proper interpretation of 
recent experiments 26,27 as well as old experiments 23 . We 
continue our discussion in Sec. VI where we derive full 
crossover functions for realistic setups covering the full 
range of chemical potentials with a particular emphasis 
on screening effects. We close with the conclusions in 
Sec. VII. We have relegated a technical discussion of the 
full scattering kernel as well as the matrix elements of 
the scattering matrices to the Appendix. 



II. THE MODEL, TIME SCALES, AND 
SCREENING 



A. The model 

The model Hamiltonian consists of two copies of 
the free graphene Hamiltonian for the active and pas- 
sive layer, respectively and interactions within and in- 
between layers. It reads 



(3) 



a/p 



where a denotes the active layer and p the passive. H 

denotes the free Hamiltonian in both layers, H^J t p the 
interaction within each layer, while H°^ t describes the in- 
teraction between layers. Disorder is implemented within 
each layer via -ff;„ t . The non-interacting Hamiltonian 
reads 



N . 

■Y, M 2 *[*/ f V -/*)*}] , (4) 

f=i J 



with the Fermi velocity vp, f = 1,...,N counting the 
flavors, and being the chemical potential of the indi- 
vidual layers. In the case of graphene we have N = 4 
due to valley and spin degeneracy, while for a topological 
insulator we would rather expect N — 1 (or more gen- 
erally an odd integer). The spinor representation of the 
wave-function has the following Fourier decomposition 



(x,t) = J 



d 2 k 
(2^) 2 



-if 

"2 / 



(M) 
(k,t) 



ik • x 



(5) 



where the operators c\/ 2 f arc the electron annihilation 
operators on the two different sublattices for flavor index 
/ and in layer i = a, p. We note that in topological insu- 
lators the spinorial components do not refer to the sublat- 
tice but rather to the spin degree of freedom accounting 
for their helical nature. The formulation of transport is 
simplest in a basis which diagonalizes the Hamiltonian 
Hq. This is accomplished by a unitary transformation 
from the Fourier mode operators (c\j,c 2 ^) to the basis 



of electrons and holes ( 7+0)7 _ ): 

4/(fc) = ^ fc ( 7 ; / (k)-7i / (k)). 



=( 7 ; / (k)+ 7 i / (k)), 



(6) 



Above, we introduced the complex number K by the re- 
lation 



K=k x +ik y , where k = (k x , k y ), 



(7) 



and k = |k| = \K\. Expressing the Hamiltonian Hq in 
terms of we obtain 



N r H 2 k 



(8) 



The distribution functions of electrons and holes (±) 
in the layers i = a/p read 



^(k,t) = ( 7 l t / (k,i)7l / (k,t) 



(9) 



There is no sum over / on the right hand side, and we 
assume the distribution functions to be the same for all 
valleys and spins, which is why we drop the index / from 
now on. In equilibrium, i.e., in the absence of external 
perturbations, the distribution functions are Fermi-Dirac 
functions 



zi(k,t) = f° x ( VF k) = 



i 



i 



(10) 



The current can be expressed in terms of the electron- 
and hole-operators and decomposes into 



(11) 



with 



and 



J ^ = ^E/(0^T 7 ^ (k) ^ (k) ' (12) 



Jjj = —ievp 



d 2 k (z x k) 



(2tt) 2 k 
7 t. (k)7_ (k)-7l (k)7 +0 (k) 



(13) 



where z is a unit vector orthogonal to the x, y plane. 3j 
measures the current carried by motion of the quasipar- 
ticles and quasiholes — notice the A prcfactor, indicating 
that these excitations have opposite charges. The oper- 
ator Jjj creates a quasiparticlc-quasiholc pair (it corre- 
sponds to the so-called Zittcrbcwegung, see Ref. 1 ) and is 
the part which determines the optical conductivity. For 
the purpose of this paper we can neglect its influence on 
transport properties, since we are interested in d.c. trans- 
port properties. In the framework of the Kubo formula, 
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which fully accounts for the off-diagonal parts, it was 
shown that this leads to numerically identical results 40 . 

In a particle-hole symmetric situation a current car- 
rying state with holes and electrons moving in opposite 
directions has a vanishing total momentum, and the cur- 
rent can decay by creation or annihilation of particle hole 
pairs, without violation of momentum conservation. This 
is the physical reason why at the particle hole symmetric 
point, i.e., at vanishing deviation of the chemical poten- 
tial from the Dirac point, the d.c. conductivity is finite 
even in the absence of momentum relaxing impurities. 
However, as we will see below, at finite deviation from 
particle hole symmetry a driving electric field always ex- 
cites the system into a state with finite momentum which 
cannot decay. This entails an infinite d.c. conductivity 
(even though drag can be finite), and consequently im- 
purities have to be taken into account. 



B. Sources of current relaxation, time-scales, and 
screening 

In the following we discuss three important ingredients 
for our subsequent discussions, which are disorder effects, 
interaction effects, as well as screening properties in two 
dimensional Dirac systems. 



1 . Sources of current relaxation 

Within this work we study the interplay of three 
different sources of current relaxation, which are in- 
tralayer Coulomb interaction, interlayer Coulomb inter- 
action, and disorder. 

The 1/r intralayer Coulomb interaction assumes the 
form 



N 



H int 



>— ^ f d 2 ki d 2 k2 d 2 q 

fh^kJ www 

1 A2 A3A4 

(ki,k 2 ,q)7l f 4/ ,(ki +q)7* t 3/ (k 2 
x 7i 2 /(k 2 )7i 1 /'(k 1 ). 
Here the scattering matrix element reads 



(14) 



q) 



JA 1 A 2 A 3 A 4 l k l 1 k 2,qj - 5 X 



I + A1A4 



{K{+Q*)K X 



1 + A 2 A 3 



(15) 
{K* 2 -Q*)K 2 



|ki +q|fci 

where w kl , q = w F (A 4 |ki + q| - Ai|ki|), and 

27re 2 



q|&2 



V(q,u>) 



e r q 



(16) 



is the dynamically screened Coulomb interaction. In this 
expression e r is the dielectric constant of the adjacent 
media. Note that we have neglected the scattering be- 
tween valleys since it connects points in the Brillouin 



zone which involve large momentum transfers and conse- 
quently are strongly suppressed. The two layers are at a 
vertical distance d (in z-direction) and consequently the 
unscreened interlayer Coulomb interaction reads 



U(r) oc 



1 

Vr 2 + d 2 



which after Fourier transform assumes the form 



I7(q,w) 



27re 2 

£r|q| 



-qd 



(17) 



(18) 



The Hamiltonian H-^ t which connects the two layers as- 
sumes the following form in the basis of electrons and 
holes 



N 



H ap 

"int 



E 

/,/' = ! A1A2A3A4 



E 



d 2 ki d 2 k 2 d 2 q 



(2tt) 2 (2^) 2 (2^) 2 

x T\ 1 \ 2 \ 3 \ 4 (k! , k 2 , q) 7 ^ /( (k! + qhH f (k 2 - q) 
X7L / (k 2 )7A a 1 /'(ki) 



(19) 



with 



JAiA2A 3 A4l k l) k 2,q) - 5 X 



1 + A1A4 



(gj + Q*)gi 
|ki + q\h 



1 + A 2 A 3 



(20) 

{K* 2 -Q*)K 2 
|k 2 - q|fe 2 



In order to discuss situations away from the Dirac point 
we have to include the effect of disorder, which is required 
in order to obtain finite individual layer conductivities. 
This is required since at finite chemical potential an elec- 
tric field excites a finite momentum state, which can only 
be relaxed due to translational invariance breaking. We 
assume the following form of the disorder potential 

hl - E / ^(x)*^)*; (x) , (21) 



with 



v dis (x) = ]T^ 

t-r 1 ex — x 



(22) 



Here Xj denotes the random positions of charged impuri- 
ties, assumed to be close to the graphene sheet, having a 
charge Ze and average spatial density Pi mp . The disorder 
Hamiltonian i?di s in terms of the 7^ reads 



i a—1 X 1 A; 



x exp[bq • (ki - k 2 )]7* t i/ (k 1 )7^ / (k 2 ), 



where 
U Xl x 2 (k 1 ,k 2 ) = 



2nZe 



e r |ki 



__ I 
W\ 2 



1 + AiA 2 



K{K 2 
k x k 2 



,(24) 
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which corresponds to unscreened Coulomb scatterers. 
Note that even though we compute specific results for 
Coulomb interacting particles and Coulomb impurities, 
the formalism easily generalizes to arbitrary isotropic two 
body interactions and disorder potentials coupling to the 
local charge density. In the following we assume that dis- 
order only acts within one layer and remains unscreened 
even at finite chemical potential. This does not influence 
any conclusions drawn from our analysis and using scalar 
impurity potentials would yield identical results. 



2. Time scales 

The transport timescales within a layer have been dis- 
cussed before 10 and we repeat the major results here. For 
a clean system at the Dirac point we find that electron- 
electron interactions induce a finite inelastic scattering 
rate. Introducing the 'fine structure constant' 



a 



e r v F 



(25) 



which has a logarithmic scaling 4 we find that close to 
zero doping it is on the order of 



-1 2^bT 



h ' 



(26) 



and thus essentially set by the temperature. This is 
a hallmark of the quantum criticality of the undoped 
graphene system 6 ' 8 ' 9 . The full crossover from quantum 



critical to Fermi liquid is described by 

, 2 k B T 2 /h 
max[/c bT, fj] 



~ a 



(27) 



where at larger doping, when the chemical potential fj, 
exceeds fcsT, the inelastic scattering rate tends to the 
expected Fermi liquid form r" 1 <~ T 2 / 'fj,, if screening is 
taken into account. The elastic scattering rate due to 
static charged impurities is naturally proportional to the 
density of impurities and in general reads 



imp 



1 (Ze 2 /e r ) 2 p imp 
h max^T, /i] 



(28) 



We note that the inelastic scattering rate decreases with 
temperature, while the elastic scattering rate increases. 
The latter is due to the fact that low energy particles are 
more intensely scattered by Coulomb impurities. Again, 
it is worthwhile mentioning that the physics of electron- 
hole puddles is beyond this description and our results 
do not apply in the inhomogeneous regime. 



3. The effect of screening 

We introduce the two independent polarization func- 
tions II a and lip for the active layer and the passive layer, 
respectively. The random phase approximation (RPA) in 
the basis of intra- and interlayer interactions leads to the 
following Dyson equation 18 



V aa (q,u) U ap (q,uj)\ = (V U 
U ap (q, w) V pp (q, w) J \UV 



V u 
U V 



which can be solved in an elementary way yielding 



Vaa U a 



-n a (q,w) 

-Hp(q,w) 



(i + vn a ) (i + vii p ) - u 2 n a ii p 



V aa (q,uj) U ap (q,uj) 
U ap (q, w) Vpp (q, ui) 



V + (v 2 - U 2 ) Up u 

u v + (v 2 — u 2 ) n a 



(29) 



(30) 



\n\/T <C 1: The non- degenerate limit. 
A peculiarity of a theory of massless Dirac fermions at 
the charge neutrality point is the absence of standard 
Thomas-Fermi screening. This can be rationalized from 
the absence of density of states at the Fermi level. This 
means only a thermal density of states enters. The zero 
temperature polarization in the Matsubara formulation 
reads 1 ' 41 - 42 



n , P (q,w n ) 



Nq 2 



l^v 2 F q 2 +u 2 n 



(31) 



From the static limit 7r(q, ui n — 0) oc |q| the absence of 
screening immediately follows. Taking into account the 



thermal density of electrons the polarization reads 

^ / „ nn iVmax[T,d AHq| 
n a>p (q, w„ = 0, T, Map ) w [ 9 ' /J + ' 



2kv 2 f 



NT N\q\ 



2'KV 2 p l&Vp 



(32) 



where T plays the role of the Thomas-Fermi screening 
momentum. Since typical momenta involved in the scat- 
tering process in this regime are on the order T/vp we 
conclude that the screening only makes a small contribu- 
tion. This contribution is again controlled in a which is 
small and thus to leading order can consistently be ne- 
glected. It turns out that in the hydrodynamic regime 
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screening must only be taken seriously if one wants to 
go beyond the two-mode approximation. However, then 
the dynamic part is only important to cut off the forward 
scattering singularity 8 ' 9 . 

|/ti|/T 3> 1: The Fermi liquid regime. 
In this limit only one of the two charger carriers matters 
and consequently one can carry out a simplified analy- 
sis 18 . The real part of the retarded polarization in the 
static limit is given by the density of states at the Fermi 
level and consequently reads 

Ren a . p (q,w„ = 0,T, Ma , p ) « . (33) 

ZITVp 

This will account for the static screening in the Fermi 
liquid regime and provides the standard Thomas-Fermi 
expression for the screening vector gxF- We will use this 
approximate form in the discussion of the Fermi liquid 
regime where we use it in U ap in Eq. (76). 

In order to understand the Fermi liquid regime and 
its limiting behavior starting from the analytical formula 
Eq. (83) we also need the imaginary part of the retarded 
polarization, which is given by 

Im^^^f— 0(v F q-\aj\) . (34) 
2ttv f vfQ 

We find that the correct polarization functions matter 
for both the determination of the correct distance and 
temperature behavior of the drag resistivity as well as 
for the orders of magnitude in the drag resistivity when 
compared to experiment. 

III. THE KINETIC APPROACH FOR THE 
BILAYER SYSTEM 

The Boltzmann equation approach has been used in 
the context of single-layer graphene in the collision- 
dominated hydrodynamic limit before 9 . We assume that 
the quasiparticle description remains valid in all regions 
of interest in our problem. The equation of motion for the 
quasiparticle distribution function schematically reads 

d t f - F ext dk/ = -7 coll (35) 

where / is the distribution function, F cxt denotes the 
external force, dt accounts for some temporal variation, 
and / co n is the collision integral. 

In our case, the system under investigation has the 
generic form shown in Fig. 1 and consequently requires 
to extend the formalism of the single layer to also account 
for the presence of the passive layer and interactions be- 
tween the two layers. This leads to a total of four coupled 
equations of motion which have to be solved simultane- 
ously. Again, we can restrict our analysis to only include 
the diagonal parts of the distribution matrix. This is jus- 
tified since we are only interested in d.c. properties. For 
optical properties this would not be justified. 



A. Coupled kinetic equations 

The general structure of the coupled Boltzmann equa- 
tions in the stationary limit, d t f = 0, assumes the form 



eEV k /^ = 


raa 
~ 1 C 


-Ic P 


raa 
- J dis 


eEV k /! = 


raa 

~ 2 c 


T ap 
1 G 


raa 
-Mis 


= 


TPP 

1 C 


T pa 
1 C 


tPP 
J dis 


= 


TPP 

1 C 


T pa 
1 C 


tPP 
1 dis 



The two uppermost lines account for the active layer in 
which both electrons and holes are subject to an applied 
electrical field. The lower two lines account for the pas- 
sive layer, in which no field is applied requiring the left- 
hand side to be zero. There is a number of collision 
terms, where Ig° and If-P account for the scattering due 
to Coulomb interaction within a layer (a and p stand for 
active and passive layer respectively), I^F and ac- 
count for inter-layer scattering, while 1%? and I^f s denote 
scattering due to disorder within the individual layers. 
The explicit form of the collision terms is presented in 
Appendix A while the matrix elements of the scattering 
matrix are defined in Appendix B. Scattering between 
the active and the passive layer only includes processes 
which arc of the density-density (large-TV) type. One 
could faithfully describe this by an effective plasmonic 
mode for the passive layer coupling to the active layer 
thereby reducing the number of degrees of freedom 43 . 
However, we choose to work in the basis described in 
Eq. (36). The effect of drag can easily be understood 
from the Boltzmann equation. In linear response the dis- 
tribution functions in the active layer /± are driven out 
of equilibrium linearly in the applied field. Consequently, 
we have to plugging this ansatz into the lower two lines 
the term I^P. This implies that now the lower two lines 
also include a part which is linear in the applied field. 
This indirectly serves as a 'source term' for the distribu- 
tion functions f± in the passive layer. In order to solve 
the lower two Boltzmann equations in linear response it 
follows that we now have to choose the deviation of /± 
from equilibrium to also be linear in the applied field in 
the active layer. In a Kubo formula approach this ef- 
fect is captured by the standard Aslamazov and Larkin 
diagrams 44 . 

B. Variational ansatz and choice of modes 

As discussed in Sec. Ill A, the distribution function 
of the quasiparticles in both layers have to be expanded 
to linear order in the applied electrical field and conse- 
quently assume the schematic form 

f a/p = f ± a/p (vFk) 

eV F k a,a/p ( f 0,a/p\ a/p ( V F k\ 

(37) 
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which provides the starting point of the subsequent dis- 
cussion. The solution strategy is to choose an ansatz for 
the functions x± P which is related to the slow modes in 
the problem. In the non-degenerate limit the analysis re- 
quires only one mode to yield an asymptotically exact re- 
sult 9 . In the degenerate limit with fi/T ^> 1 the most im- 
portant mode is the momentum mode. Both modes share 
the property that they can annihilate the divergence in 
the forward scattering amplitude of the Coulomb colli- 
sion kernel, which is a peculiarity of electrons with linear 
dispersion in two dimension 8,9 . These modes thus con- 
stitute the leading contribution to current relaxation to 
leading logarithmic accuracy, which has been used before 
to describe the crossover of the single layer conductivity 
from the non-degenerate limit to the degenerate limit 10 . 
For a in depth discussion of the logarithmic singularity in 
forward scattering we refer the reader to Ref. 10 . The ap- 
propriate minimal ansatz for our purposes consequently 
is given by 



a/p ( V F k\ 

xl [-jr) 



+v a /p _1_ , , a /p Vpk 
±Xo± + Xl± -jT 



(38) 



where xa is associated with particle-hole symmetry and 
the ± accounts for that while Xi refers to the momentum 
conservation. The mode xo dominates transport in the 
non-degenerate limit. On the other hand \\ dominates 
in the degenerate limit. The solution of the Boltzmann 
equation now is equivalent to determining the coefficients 
Xo± an d Xi± ') which can be mapped to a matrix inver- 
sion problem. The general formalism is an application 
of the variational principle for coupled Boltzmann equa- 
tions 39 which is explained in great detail below. 



C. Matrix formalism for drag 

Using the set of functions defined in Eq. (38) the Boltz- 
mann equation and its solution can be mapped to a 
matrix inversion, where the matrix acts in a combined 
space of layer indices, electrons, holes, and modes. This 

gives access to the expansion coefficients Xo^± an d Xi^±> 
which then allows to determine the individual and trans- 
conductivities. The major numerical effort within this 
approach goes into a faithful calculation of the matrix 
elements of the collision kernel. The space of func- 
tions, layers, and particle nature allows to define a vector 

X = (xo+>Xo->Xi+>Xi-,Xo+'Xo-'Xi+,Xi-) where the 
indices are chosen as in Eq. (38). The space of functions 
is defined by 



k k vph vpk k 



k wpk vpli 

jk|' ~T~' ~T~ 



(39) 



with i = 1,...,8. One can expand the right hand side 
collision operator in Eq. (36) to linear order in the field 



E, which leads to the following schematic expression 

-^|/!'°(i-/!' a ) = -(c aa + c ap + cz) 
o = -{c pp + c pa +c pp s ) 



Xq/i,± 



a/p 
0/1, ± 



o = -(c pp + c pa +CZ) 



a/p 

Xo/i,± 

a/p 

Xo/i,± 



(40) 



We define the scalar product between two objects in 
this space as 

(a\b) = J ^a(k)Kk) • (41) 

In the following we allow the more general case of applied 
fields in both layers. This is a generalization giving access 
to all quantities within the conductivity tensor including 
the passive layer conductivity, which in principle can be 
different from the active layer. We define a vector of the 
driving term as 

D = (3 a ,D p } with 

D a = {{e l \D a + ),{e 2 \D a _),{es\D a + ),{e i \D a _)), 

D P = ({e 5 \D p + ),{e 6 \D p _),{e 7 \D p + ),{e 8 \D p _)) (42) 



where 



Da+/p = ev^ f o, a/p ^ /r / P ) and 

D a,p = _e_V^E f J a/p^_ f _,a/p^ ^ 



Equivalently, the elements of the collision matrix are 
given by 

Cij - {e t \C aa + C ap + CZ + C pp + C pa + CZ\e 3 ) (44) 

where the superscripts on the right hand side indicate 
that the matrices act within layer space. Finally, the 
Boltzmann equation can be cast in the form 



D = C-x. 
A straightforward matrix inversion 
X = C~ 1 -D 



(45) 



(46) 



allows to determine the expansion coefficients. The pro- 
jection to obtain the respective single-layer conductivities 
and the transconductance is done via 



X a - C~ 1 -{D a , 0,0, 0,0) and 
x p = c- 1 - (0,0, 0,0,5") . 



(47) 
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The conductivities in the individual layers now read 
- a = ^E^T| |EM (48) 

and 

ffP ^E^| |EN (49) 
while the off diagonal drag conductivity reads 

- d = ^E^7 +4 | |EM , ( 5 °) 

i— 1 

where D was introduced in Eq. (42) , and N is the number 
of valley and spin degrees of freedom. We will give a 
concrete example of this formalism in a reduced setting 
in Sec. IV C and Sec. V. 



IV. NON-DEGENERATE LIMIT: \n a/p \/T «: 1 

In this limit the difference from the standard Fermi liq- 
uid behavior of Coulomb drag is expected to be largest: 
not only the transconductivity, but also the conductiv- 
ity of the individual layers are expected to possibly be 
dominated by cither inelastic or elastic scattering. This 
implies that we expect drastic changes as the ratio of dis- 
order to interactions is changed. We will find that this 
ratio can alter the drag resistivity by orders of magni- 
tude. In the discussion of this limit we neglect the effect 
of screening due to the lack of density of states, which is 
the rational given in Sec. II B 3. Our analysis does not 
capture the regime of electron- and hole-puddles and we 
assume in the following that the temperature T is high 
enough to be beyond the inhomogeneity scale. With in- 
creasing sample quality we expect that this regime can 
be pushed to rather low temperatures thereby increasing 
the domain of validity of our analysis. All calculations 
throughout this section have been performed explicitly 
for graphene, meaning N = 4, meaning all numbers are 
calculated for this case. Since we are going behind large- 
N including crossed diagrams in our analysis we can not 
deduce the result for different values of N by simply scal- 
ing. 

A. Both layers at the Dirac point 

If both layers are at the Dirac point there is no drag 
due to particle-hole symmetry in both layers: Plugging 
the parametrization Eq. (37) and Eq. (38) into the Boltz- 
mann equation of the passive layer it is straightforward 
but rather tedious to show that the equilibrium distribu- 
tion function 

mKt) = £ p (v F k) (5i) 



solves the Boltzmann equation for the passive layer. This 
immediately implies that in the setup shown in Fig. 1 we 
have V2 = R2 = 0, meaning the aforementioned absence 
of drag. However, in this limit we have a well-defined 
single-layer charge response even in the absence of disor- 
der since the current carrying state in the active layer is 
of zero total momentum. This was found to be given by 

a = 0.76^ (52) 
by a numerical solution of the Boltzmann equation in 

q 2 

the leading logarithmic approximation 9 with a = 
being the dimcnsionlcss fine structure constant. We note 
that a similar calculation taking into account only the 
large-N diagrams (discarding crossed diagrams) was car- 
ried out by Kashuba 8 and the result was a' = 0.69^2. 
In the limit Td/vp > 1 we expect and find that the ac- 
tive layer conductivity extrapolates to this isolated sin- 
gle layer conductivity, an. However, we find there is a 
substantial effect of the passive layer on the transport 
properties of the active layer in the limit Td/vp < 0(1). 

Overall, in the absence of screening effects the conduc- 
tivity of the active layer (in the passive layer it is zero) 
is a one parameter function of the type 

a a (T,d)=a a (™) . (53) 

We calculate this crossover function from a full numeri- 
cal solution of the four coupled Boltzmann equations (at 
charge neutrality one could in principle reduce the num- 
ber by a factor two since there is a particle-hole symmetry 
to exploit). 

It is shown in Fig. 2 function of the dimcnsionlcss 

parameter Td/vp. It interpolates from 0.36^2 at low 

2 

temperatures to 0.76^2 at high temperatures, which is 
the isolated single layer conductivity, see Eq. (52). This 
implies that for a fixed distance as a function of tempera- 
ture the single layer conductivity can change by roughly a 
factor of two. We can easily rationalize the results found 
numerically in the limit ™ — > 0: in the absence of charge 
currents between the layers the charge carriers in the ac- 
tive layer scatter from charge carriers in the active layer 
as well as from those in the passive layer. Both types 
of processes share the same Coulomb potential due to 
^ — > 0, which implies the exponential screening factor 
in Eq. (18) is not active for typical momenta. The scat- 
tering times associated with inelastic scattering within 
and in-between layers will thus just add up. We can guess 
the result for the conductivity from the two single layer 
results mentioned above: For a clean sheet of graphene 

q 2 

the conductivity 9 assumes the value 00 = 0.76^2. This 
result takes into account all diagrams of the Born ap- 
proximation, thus also crossed diagrams beyond large- 
N. In a calculation which discarded crossed diagrams 
and concentrated on large-N diagrams Kashuba 8 found 
the conductivity was given by CTq = 0.69^2 for a sin- 
gle layer. Since scattering across layers only involves 
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conductivity a a as a func- 

The conductivity interpolates be- 
tween two limiting values at high and low temperatures 
and is described as a universal function which only depends 
upon Td/vF- For low temperatures the behavior is given 



FIG. 2: Minimal d.c. 
tion of T in units ^F-. 



by lim 



Td/vt 



+o (a (Td/v F ) 



0.36 T 



temperatures it is given by lim 



Td/v f 



(cro 



while for large 
-a(Td/v F )) oc 



density-density type scattering (plasmons) it is faithfully 
accounted for by large-N type diagrams and we expect 
that in the limit Td/vp <C 1 the conductivity is given by 



0.36, 



ha 2 



(54) 



which corresponds to adding the inverse scattering times. 
The low temperature behavior can be rationalized from 
the exponential factor in the interlayer potential: the typ- 
ical momentum of electrons and holes involved in the elec- 
tronic transport is the thermal momentum q typ = T/vf- 
Expanding the exponential of inter-layer Coulomb in- 
teraction for — <C 1 at this typical momentum thus 



yields a correction U(q,T) 
V(q,T) (l 



V(q,T)(l-q typ d) 



Td 

IT 



) , which suggests a linear variation of the 

conductivity for very small temperatures. 

An interesting question is to speculate whether this 
effect should be visible in available geometries. Typical 
values of layer separations which are currently used in 
experiments are on the couple of nanometer range. For 
a distance of d w 10 nm this implies that the typical 
crossover temperature which we extract from demand- 
ing that ^ « 0.2 - 0.4 is given by T cross = *$■ « 
150^ — 300^. Consequently, in experiments one can eas- 
ily reach the temperature range where interaction effects 
are enhanced. 

We note that we performed the same calculation by 
describing the passive layer in terms of plasmons. It 
was shown before that due to particle-hole symmetry this 
mode remains in equilibrium and thus one can calcu- 
late the scattering of electrons and holes from equilibrium 





iu d2 (q) 
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k tUdl(q) 
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FIG. 3: Possible experimental setup which can effectively in- 
crease the role of inelastic scattering close to the Dirac point. 
The layers adjacent to the central active layer, in which a cur- 
rent I is driven serve as reservoirs which provide a source of 
inelastic scattering. 



bosons which simplifies the analysis. We have checked 
that the results are in perfect agreement. 



B. Increasing inelastic scattering close to the Dirac 

point 

The above results immediately bring about a promising 
route towards increasing the effect of inelastic scattering 
in experiments carried out in the vicinity of the Dirac 
point in ultra clean samples: stacking a larger number 
of monolayers effectively increases the effect of inelas- 
tic scattering since the inverse scattering times due to 
interactions of the individual layers add up as long as 
-C 1. Current experiments on bilayers have demon- 
strated that distances d = 1 nm are conceivable without 
leakage currents which implies that one could arrange a 
large number of layers in a sandwich structure and still be 
well below the crossover scale for temperatures up to the 
one hundred Kelvin range. A possible schematic setup 
is shown in Fig. 3 where the central layer is the active 
layer in which the current is driven while the surround- 
ing passive layers solely increase the effect of interactions 
but remain in equilibrium themselves. We propose that 
such a sandwich structure can facilitate experiments in 
the limit of the sought after hydrodynamic interaction 
dominated regime 10,12 . 



C. Finite chemical potential 

In the following we concentrate on equal chemical po- 
tentials (our analysis trivially includes the situation of 
equal doping but with different types of charge carriers 
in the layers, meaning fi a = —fi p , which results in an 
overall minus sign) and comment on different chemical 
potentials in the individual layers in Sec. VI. Equal dop- 
ing can be achieved in samples in which the charge carrier 
density can be controlled individually by separate gates. 
As discussed above, when both layers are at zero doping 
there is no drag. Even though for finite doping this is 
not true any more, the electron and hole density in this 
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regime still is mainly thermal. Like in the case discussed 
in Sec. IV A we expect to find a strong interaction effect 
and enhanced inelastic scattering here as well, especially 
in the limit — <C 1, where the intcrlayer interaction is 
essentially undamped. 

The conclusions of the following discussions hold for 
arbitrary finite chemical potential but for small chemical 



potential we can carry out a simplified analysis. In this 
regime we can neglect the effect of screening and one can 
relate different parameters in terms of disorder and in- 
teraction strength easily by scaling the matrix elements. 
For simplicity, we also assume that in both layers there 
is the same amount of disorder. The collision matrix has 
the aforementioned block structure and assumes the form 



a 



C aa (Ti)+C a P(-p,d) + ^C dis (ji) 
-C ap (]i,d) 



C aa (jI)+C a P(Ji,d) + ^C dis (jl) 



= a 2 C,,. 



(55) 



where a = ^ Ze 2°/ t y p _ — , d — , and \i = !p are dimen- 
sionlcss parameters and i, j denotes the layer indices a 
and p. The dimensionless parameter a 2 corresponds to 
the ratio of elastic scattering to inelastic scattering, i.e., 
Timp/Tee- Furthermore, we have the driving terms 



D a (p) = D p (p). 



(56) 



From this we find rather simple expression for the indi- 
vidual conductivities 

cr a (a,a,d,~p) = cr p (a,a,d,Ji) 
Nire 2 



ha 2 
Nire 2 



£>„(/!) -C" 1 ■£>„(# and 



a d (ot,a,d,n) = -^-D a (fi) ■ C a * ■ D a {p) (57) 



leading to 



(3 a (?) • Caa 1 ■ D a fc)) - (3 a (/Z) ■ Cap 1 ■ D a fc)) 

o h 



Nne 2 



g(a,d,n) . 



(58) 



From the above collision matrix Eq. (55) it is obvious 
that we have two limiting cases: a = corresponds to 
the disorder dominated limit, while a — > oo corresponds 
to the clean limit. 

The disorder dominated limit: a — > 
Using the collision matrix and its components introduced 
in Eq. (55) it is straightforward to find that the leading 
order in a expression of drag reads 



p d ~ a 



h D a (jl) ■C d ; i l(jl)C ap (d,Ji)C^(ji) ■ D a {-p) 



Nire 2 



(59) 



This expression is well behaved and no singular matrix 
operations are involved. This is due to the fact that in the 
disordered limit the presence of impurities breaks trans- 
lational invariance and the individual conductivities are 



always well defined. The above expression is consistent 
with the standard approximation 



Pd = 



-o- d 



>a>Jp 



0~aO~p 



(60) 



This approximation does not hold in the interaction 
dominated limit which in contrast to traditional two 
dimensional electronic systems might be attainable in 
graphenc. The interaction dominated limit: a — > oo 
For zero chemical potential the conductivity in all lay- 
ers is well defined even in the absence of disorder due to 
particle-hole symmetry. However, drag also vanishes for 
the very same reason. This changes for finite chemical 
potential: in the limit of vanishing disorder the individ- 
ual layer conductivities as well as the transconductivity 
diverge. However, it turns out that the drag resistance 
can still be finite. This is an effect of the boundary con- 
ditions which are such that no current is allowed to flow 
in the passive layer. The full collision matrix was intro- 
duced in Eq. (55) and we see that in the limit of vanishing 
disorder (a — > oo) it assumes a simplified form. At finite 
chemical potential it turns out that the matrix C aa {jt) is 
not invertible due to the existence of the momentum zero 
modes, which are excited. This is not true for C ap (~fi, d) 
which is invertible. However, if C aa (~p) = the full matrix 
is not invertible. This implies we have to take care per- 
forming the limit a — > oo and should not do so from the 
outset. Consequently, we need the effect of disorder act- 
ing on the zero modes of C aa (]l) in order to regularize the 
response within the individual layers. In order to isolate 
the space of zero modes we first transform the collision 
matrix by the matrix U which diagonalizes the collision 
matrix C aa (p), meaning we perform the operation 



& 



U 
u 



u- 







u- 1 



where U is chosen such that 



(61) 



(62) 



where V aa (p) is a diagonal matrix with zero eigenvalues 
corresponding to momentum conservation. In the follow- 
ing we will make the simplifying assumption (only for the 
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purpose of concise presentation) that the sectors of the 
zero modes and the other modes do not mix. This implies 
we can carry out a simplified discussion of the individual 
conductivities. We introduce an index which refers to 



the the space of zero modes, and 1 referring to the other 
components. The transformed collision matrix assumes 
the form (remember 2?qq = 0) 



C = a 2 





V o 







-C$(fi,d) 







\ 



p^) + Ci7GM) + ^Ci d rG«) o 

o c a p GM) + ^ s (m) o_ 

-Cll{H, d) VfiQi) + C?f (m, d) + ±C*?(ji) J 



Using this expression one can first analyze the individual conductivities which read 



Nire' 
ha 2 

Nire 2 
ha 2 



D a (jl)U 



-i 2 



(c$r) 1 o 



1 UD a (-p) 



(C s ) 1 o 



4(p) 



(63) 



(64) 



and 



./We 2 / it),, 71 ,.\ / *2 (^oo s 



^ aW, awn Q (^ii + cii-ciiivtt + ctiy 1 ^ ciiivit + cii)' 1 




(65) 



We can now understand why the conductivities diverge. In the case of interlayer conductivity a a there now is one 
contribution cx dP a (ji)^ (C$ s ) _1 which diverges in the limit a — > 00, as it should. The very same term is 

responsible for the divergence of Od- Despite these divergencies the transresistivity remains finite, since the most 
severe divergences cancel. We have done this explicitly and found that 

2 2h 1 

Jim pd = -a — — ^ : r=ri : , (66) 

Nne di(-p) (v% + cn en (Vtf + cur 1 eg) (1 - c\i (pit + erf)" 1 ) 4(7*) 



which indeed is finite in the limit a — > 00 since no 
singular matrices are involved and all information about 
disorder is gone. However, the validity of the above ex- 
pression is restricted to finite values of d, which becomes 
apparent from the fact that if C ap — > in the limit d — > 00 
it remains finite. This is rooted in the implicit assump- 
tion of the above analysis that Cqq 3> =zC dls , which does 
not hold in the limit d — > 00. The above finiteness of 
the response is similar to the so-called universal conduc- 
tivity, where the Kubo expression is regularized by finite 
disorder, which drops out in final expression, allowing to 
extrapolate to the clean limit. 

Overall, we have shown that drag resistivity is indeed 
described by a function of the type Eq. (58) which is 



always finite 

0< \g(a = 00, d < 00, /x)| <oo, (67) 

meaning the drag resistivity remains finite even in a 
clean system (strictly speaking we used /Z <C 1 to dis- 
card screening without loss of generality). For very small 
values of fi/T we expect the drag to be cx {p-/T) 2 for 
symmetry reason which is also backed up by our nu- 
merical analysis. We studied drag in the vicinity of the 
Dirac point as a function of the dimensionless parame- 
ter a in great detail. The clean system is found in the 
limit a — > 00 . A first observation is that for finite chemi- 
cal potential we can make the extrapolation to the clean 
system and find that the drag resistivity remains finite. 
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This can be seen in Fig. 4 where for distances d = 0, 1 
and for fi/T — 1/20 we have plotted the function 



?(a,d,/i) = l/p 2 g(a,d,p.) 



(68) 



as a function of a. We have checked that this extrapo- 
lation can be performed for any finite chemical potential 
and the limiting value in the clean system, a — > oo, in- 
creases upon decreasing the chemical potential. 




FIG. 4: Crossover function describing drag all the way from 
the disorder limited (a — > 0) to the clean system a — > oo 
for different dimensionless distances d — 0,1 and Ji = 1/20. 
Importantly this quantity saturates which we have checked 
explicitly for different realizations of ~fi. 




FIG. 5: Drag as a function of p/T for different disorder 
strengths ranging from the disorder limited (a — > 0) all the 
way to the clean system a — >■ oo. We find that in the limit 
p — y we end a with a finite drag upon extrapolation. The 
value to which pd extrapolates is given by pd — — ^ where ao 
is the single layer conductivity and was defined in Eq. (52). 



fitting curve 

A* 2 

- Pd /{a 2 h/e 2 ) = — i ■ (69) 

ai(a ) + a 2 (a )p + a 3 (a )pr 

The results of this fitting procedure for ai(a 2 ), a 2 (a 2 ), 
a^(a 2 ) and are shown in Fig. 6. 



This is in contrast to the statement based on particle- 
hole symmetry that drag at the Dirac point is zero. This 
signals that in the limit a — > oo, which corresponds to 
the ballistic system, extrapolating to zero density must 
become singular which we show explicitly. We have fur- 
thermore studied the drag resistivity as a function of fi/T 
for a set of different disorder realizations, meaning differ- 
ent values of the parameter a. For finite disorder the 
overall shape is such that there is a maximum of drag 
resistivity. In the limit of very low chemical potentials 
compared to the maximum position we indeed find be- 
havior of the (p/T) 2 type which is consistent with expec- 
tations based on symmetry considerations. We find that 
upon decreasing disorder the maximum of drag shifts to- 
wards lower chemical potentials and pushes to zero in the 
clean limit. Consequently, with decreasing disorder the 
quadratic regime becomes increasingly small and it van- 
ishes in the limit of zero disorder. These behaviors are 
extracted from our numerical results which are summa- 
rized in Fig. 5 where again we have discarded the role of 
screening and all the curves are plotted for d = 0. 

The extrapolation to a — oo, i.e., the clean case re- 
quires some care since it is very sensitive to small numer- 
ical errors. We show below how we obtained the limiting 
curve in Fig. 5. We extract this behavior from fitting the 
numerical curves for different disorder realizations and 
extrapolating to the clean limit. In order to do so we 
have fitted the ascent of the curves with the following 
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FIG. 6: Fitting parameters ai — as as a function of the pa- 
rameter 1 /a 2 . Zero corresponds to the clean limit while large 
values correspond to the dirty limit. Most importantly, to 
within numerical accuracy we find that ai — > as a — > oo. 
This implies that in principle in the clean system at charge 
neutrality there can be finite drag. 

Most importantly, we find that in the clean limit, a\ — > 
within our numerical accuracy, while a 2 and 0,3 remain 
finite. In the limits of /i/T <C 1 and 1/a 2 C 1 we find 

( 1 \ a 2 h JI 2 

p d /K 1, < 1 = ^ ^- r . (70) 

H V " / e 2 0.76/Z 2 + 0.012^ V ; 

This implies that upon performing the limit — >■ 
after performing the extrapolation to the clean limit, a — ¥ 
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oo, we end up with finite drag. Consequently, the drag 
resistivity depends upon the order of limits according to 

lim Jim p^ = 

Jim lim p^ = , (71) 

a— >oo fj, a — tip— >0 

where <jq was introduced in Eq. (52) and denotes the sin- 
gle layer conductivity. A finite limiting value has already 
been observed in the recent work by Schiitt et al. 3e . Inter- 
estingly, studying their numbers shows that their result is 
consistent with the statement lim Ma=Mp ^ nm 5->cx> Pd = 
— even though there seem to be numerical discrep- 
ancies. Schiitt et at worked in the framework of the 
large-A approximation, in which crossed diagrams are 
neglected as opposed to our analysis. Within that ap- 
proximation the single layer conductivity has to be re- 
placed by ctq = 0.69^j 8 which is consistent with the 
numerical value found. We thus conclude by saying that 
our results are not only qualitatively but also quantita- 
tively compatible taking into account the slightly differ- 
ent approximation schemes used. 

It is important to point out that in the physical system 
where a is finite the drag at the Dirac point is always zero 
for symmetry reasons. We note that there still is a way 
towards finite drag at the Dirac point which is rooted 
in including a 3 processes which however is beyond our 
scope 36 ' 46 . In that case the symmetry arguments ensur- 
ing zero drag are invalid and finding finite drag at the 
Dirac point is possible. To summarize, our result sug- 



gests that for extremely clean samples in the ballistic 
limit it the regime of doping in which the drag resistiv- 
ity goes to zero can in principle becoming very narrow. 
However, we stress that we do not think that this effect 
is at the heart of the experimentally observed zero-bias 
drag in graphene 28 . For any finite disorder level the fact 
that the drag resistivity drops to zero at the Dirac point 
seems inevitable to the order we consider here. More 
likely it is rooted in the a 3 -contribution 36 or related to 
a mechanism which relies on energy transfer between the 
two inhomogeneous layers 37 . 



V. FERMI LIQUID REGIME: p/T > 1 

In the Fermi liquid regime screening effects become 
crucial, see Sec. II B 3. In the limit of strong doping it 
is reasonable to consider in both layers only one species 
of charge carrier. We assume that the chemical potential 
in both layers is large and positive implying that only 
electrons are involved in the processes. Since we are in 
the disorder dominated regime we restrict our analysis 
to the momentum mode implying that instead of work- 
ing with matrices of dimension eight we can work with 
matrices of dimension two. We have checked that this 
reduction in the Fermi liquid regime leads to numerically 
identical results with the calculation involving all 64 ma- 
trix elements. Again, both layers are characterized by 
identical charge carrier concentration as well as disorder 
level. The reduced Boltzmann equation reads 







(e 3 |C° 



■C a >+CS?\e 3 ) 



(e 7 \C pa \e 3 ) 



(e 3 \C° 



e 7 



(e 7 \CPP +CP a +CZ\e 7 ) 



P 

Xi,+ 



(72) 



For reasons of momentum conservation there is no re- 
laxation of the current due to intralayer interactions, 
which implies that (e 3 |C Qa |e 3 ) = (e 7 \C pp \e 7 ) = 0. Defin- 
ing (e 3 |C ap |e 7 ) = -Ccoui this implies that (e 7 \C pa \e 3 ) = 
-Ccoul, while (e 3 |C ap |e 3 ) = (e 7 \C pa \e 7 ) = C Coul for 
symmetry reasons, consistent with an infinite response 
in absence of impurities. Furthermore, we choose 
(e 3 |Qis|e 3 ) = (e 7 |Cj?|e 7 ) = C d is and introduce the short- 
hand V — (e 3 \D+). The individual conductivities read 



a v = 



°d = 



Nne 2 


f 2 (Cdis + Ccoul) 


hT 


^dis + ^ CcoulCdis 


Mire 2 


2? 2 (Cdis + Ccoul) 


hT 


^dis + 2 CcoulCdis 


Nire 2 


f 2 Ccoul 



and 



hT 



C 2 



2 CcouiCdis 



(73) 



From this the drag resistivity obtains as 



Pd = - 



hT Ccoui 
A^rre 2 V 2 



(74) 



We stress that no further approximation has been used 
to arrive at this final result. The driving term is given 

by 



V 



2ttT 2 



(75) 



15 



and Ccoui assumes the relatively simple form 
d 2 k d 2 k x d 2 q 



Ccoul 



ANv% 

~T^~ J (2tt) 2 (2tt) 2 (2tt) 2 
x 2nS (vpk + vpki — vp\k + q| — v F \k\ — q|) x 
X q-q|f ++++ (k,k 1 ,q)| 2 X 

X /° (k)f + (fciXl - f + (|k + q|))(l - /°(jk! - q|)) , 

(76) 



where T ++++ was defined in Eq. (20). This leads to the 
following expression for the drag resistivity 



, /m ^ ^ (TY2lT 2 V% 



x |f/ ap (q,w F (fc- |k + q 



(2 2 fc d 2 fci d 2 g 



(2tt) 2 (2tt) 2 (2tt) 

1+ (A- + Q)** 



2 <5 (u^fc + «f&i - v F \k + q| - v F \ki - q|) q- q 



fc|k + q| 

x /° (fc)/° ( fcl )(l - /° (|k + q|))(l - f + (|k! - q|)) 



1 + 



(gi - Q)*ffi 
fci|ki - q| 



(77) 



r 



We continue to show how our approach recovers the 
standard formulae of drag 18,19 in a Fermi liquid in a 
straightforward manner, which is usually derived in the 
framework of a Kubo formula calculation using Ward 
identities. There are two key rewritings that we use in 
the following. We use the identity 



Both expressions have to be used twice for k and k\ sep- 
arately. In this expression we use 



8{v F k + v F k\ — v F \k + q\ — v F \k\ — q|) = 
J dujS(uj — vpk + vp\k + q\)S(uj + vpki — vp\ki — q|) 

(78) 

as well the identity between Bose and Fermi functions 

fl{vpk){\- f + (vpk-L0)) = 

n B (cu) (f° + (vpk) - f + (v F k - to)) . (79) 



n B (uj) 



(80) 



which is the standard Bose function. The imaginary part 
of the retarded polarization function can be written as 



Imn++(q,w) = Nn j -^f 5 (u-v F k + v F \k + q\) 1 - ( 



(K + Q)*K 
fc|k + q| 



(/> F fc)-/°Mk + q|))) (81) 



where the superscript ++ signals we only consider the we obtain the well-known formula 18,19 ' 34,35 

electronic part (or for the hole part). This becomes 

asymptotically exact in the limit ji/T — > oo which is 
what we concentrate on. Plugging these expressions into 
Eq. (77) and using 



n B (uj)n B (-uj) 



4 sinh 2 ^ 



(82) 
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ft,(/i/T»l) = 



h 2 



e 2 7T 2 HaUpVaVp T 



7^ I du! 



d 2 q 



(2tt) 2 sinh 2 ^ 



|[/ ap (q,c))| 2 Imn++(q, W )Imn++(-q,- W ) , (83) 



r 



where v a , P is the density of states at the Fermi level, 
which in graphene is given by v ap — N J ia ^ p . The limiting 

values of the transresistivity can now be estimated upon 
knowledge of the screened interaction potential which ac- 
cording to Eq. (30) is given by 



U ap (q,uj)) 



U 



(i + vii a )(i + vu p )-um a u p 



(84) 



In the ballistic limit the polarization function reads 

Imn++(< Z , W , M /T» 1) = v a , v — Q(v F q- M) . (85) 

v F q 

Furthermore, we use Eq. (33) to describe the effect of 
screening. 

Rescaling variables to q — > ^ and lo — > ujT yields (we 

use d = ^f) and performing a number of manipulations 
we can rewrite the drag resistivity as 



Pd 



T V 



16tt 1 



(86) 



where we used kp 

aNkpd 



K 



— and introduced the shorthand 

VF 

. Furthermore, we introduced the function 



g(d, k) = 



dqq 3 e 



-2q 



o ( ( ^ + i) 2 _ e -2 9 ) J-q/d^h 



q/d dcouo 2 



.(87) 



As mentioned before, we are considering the limit \xjT » 
1, but instead we leave d as well as kpd arbitrary. This 
implies there are different regimes, which one can access 
in the above formula. We have checked that all results 
from Eq. (83) are numerically identical to the ones ob- 
tained from directly integrating Eq. (77). 



A. Drag in the limit d <C 1 

This is the limit in which we expect to recover the stan- 
dard results of Fermi liquid theory. We can approximate 
the function 



g{d<. 1 



K) « / 
JO 



dqq 3 e 



3 p -2q 



8?r 



dqq 3 e- 2q 



duuj 
sinh | 



3 ^ {[q K +lf -e- 2 *)' 



implying it is independent of T and thus drag has the 
standard T 2 -behavior. This expression has two limiting 
behaviors as a function kpd. 



k F d < 1: 

In this limit, we find that the role of k cannot be ne- 
glected and consequently we have 



g(d<€. 1,k > 1) 



87^2. 
3 ji/ k 

8?r 2 In k 

This implies the drag resistivity reads 



dqe 



-2q 



Pd 



h (T 



e 2 \/i 



8a 2 , aNkpd 

In 

3?r 2tt 



(89) 



(90) 



kpd > 1: 

In this limit, we find that the role of n can be neglected 
and consequently we have 



g{d<€. 1,k> 1) 



87T 2 r°° dqq 3 e- 2 i 
3 Jo (l-e- 2 iy 
tt 2 C(3) , 



where ((x) is the Riemann function and £(3) 
This implies the drag resistivity reads 



Pd 



h (T\ 2 16tt 3 C(3) 1 

iV 4 (kpdf a 2 



e 2 \P 



(91) 
1.202. 

(92) 



We have checked both statements against numerically 
integrating Eq. (83) and the results are shown in Fig. 7. 



B. Drag in limit d 1 

This limit only allows for one limit of k, since by con- 
struction the parameter which enters is given by dJL, 
which by construction is large in this limit (remember 
n/T >• 1). Consequently we have k«1 and taking into 
account the fact that the integral over momentum q is 
cut off on the scale one due to the exponential factors we 
find 



f c 

g(d > l,/c< 1) w / 
Jo 



,3 e -2q 



dqq 3 e 
(1 - e- 2 i) 2 



q/d 
q/d 



dw 



IT* 

15d 



For the drag resistivity this implies 

h T 16tt 5 1 
Pd ~ eVl5A^ (kpd) 5 a 2 ' 



(93) 



(94) 
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FIG. 7: Crossover from hid- to l/d 4 -behavior in the d< 1. 
The curve was obtained for kp = 10 and T = 0.01. The 
crossover takes place roughly at kpd ps 1. 



-Pd[a.u/ 




FIG. 8: Crossover from T 2 to T linear behavior. The curve 
was obtained for kp = 100 and d = 1. The crossover takes 
place roughly at Td ~ 0.2. This scale is achievable in cur- 
rently available samples. 



meaning it is linearly proportional to T and inversely pro- 
portional to d 5 . Again, we have numerically integrated 
Eq. (83) and verified these predictions, as shown in Fig. 8. 
We have found that the deviations from T 2 -behavior be- 
come visible on the scale d ps 0.1 — 0.2. Most importantly, 
this shows that the use of Eq. (85) as approximate form 
of the imaginary part of the polarization function is fully 
justified. 



C. Connection to experiments 

In the limit fi/T > 1 we have identified a variety of 
regimes depending on whether dT/vF and k F d small or 



large 



Pd 



h ( T\ 2 8oi 1„ aNkpd 
"i 2 \ u y 3?r 2-7T 

(t\ 16tt 3 C(3) 1 

\fij N 4 (k F d) 4 a 2 

h T 16tt 5 1 
e 2 ft 15N 4 (k F d) s a 2 



d < 1 k F d < 1 

d < 1 k F d > 1 ( 95 ) 
d > 1 kpd > 1 



We have studied Eq. (86) numerically, especially in the 
limit kpd 3> 1, to determine the crossover. We find that 
in this limit the deviations from the T 2 -behavior become 
sizeable already at d ps 0.1 — 0.2, leading to a broad 
crossover region. For a sample with interlayer distance 
d ps 10 nm this translates to a crossover temperature 
on the order T cross ps 150A, meaning that the deviation 
from the standard T 2 behavior should be observable in 
current samples. This is particularly interesting in light 
of the results in the recent experiment by Kim et al. 26 . 
In this work a substantial deviation from the T 2 behavior 
was found in the temperature range above 150 — 200A". 
However, we note that Kim et al. extract their tempera- 
ture behavior from the maximum of the drag while here 
we did so deep within the Fermi liquid regime. We have 
solved Eq. (77) as well as the full problem with all parti- 
cle sorts and modes also in the regime of the maximum 
of drag and found results fully compatible with the above 
discussion. We thus conclude this section by stating that 
a possible explanation of the experimental finding of a de- 
viation from the standard T 2 Fermi liquid behavior could 
be that the experiment enters the very broad crossover 
regime where the behavior crosses over to the linear in T 
behavior. It is also interesting to note that very similar 
behavior was found in one of the earliest experiments on 
two-dimensional electron gases by Gramila et al. 23 where 
for higher temperatures large deviations from T 2 were ob- 
served. We have checked that in their work the crossover 
scale is also roughly given by d = 0.2. In a more recent 
experiment by Gorbatchev et al. 28 the authors found T 2 
behavior of the drag resistivity, consistent with the stan- 
dard Fermi liquid predictions. However, instead of l/d A 
or l/d° the authors find a 1/d 2 behavior. One possible 
explanation is that the measurement takes place for val- 
ues of kpd fa 0(1), where there is a crossover from l/d A 
to l/d° behavior, see Eq. (95). Coincidentally, the in- 
termediate range might appear as 1/d 2 . An alternative 
explanation follows the paper by Kamenev et al. 18 which 
discussed the Fermi liquid regime. In the diffusive limit 
it was found that there should be a behavior which is T 2 , 
but of the 1/d 2 type. 



VI. FULL CROSSOVERS 

So far we have discussed the limiting cases deep within 
the non-degenerate limit, fi/T <C 1, and in the opposite 
Fermi liquid regime where n/T 3> 1. We conclude our 
discussions by commenting on the drag resistivity as it 
obtains in an experiment where both layers can be gated 
individually. This is relevant in experiments, where the 
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charge density in the most general case varies in the two 
layers 26 ' 27 even though more recent experiment achieve 
equal carrier densities with very high precision 28 . This 
section mainly serves to complete the overall picture and 
for comparison to experimental findings. In this discus- 
sion we consider the fully screened interaction. In the 
case where both active and passive layers are kept at 
identical chemical potential sweeping the carrier density 
results in a curve which extrapolates between the (p/T) 2 - 
behavior close to charge neutrality to the standard Fermi 
liquid regime, Fig. 9. For this plot we have kept both 
chemical potentials identical and show two different dis- 
tances d. We observe, as expected, a strong overall de- 
pendence of the order of magnitude as well as the position 
of the maximum with varying distance. Our numerical 
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FIG. 9: Numerical evaluation of the drag resistivity for equal 
chemical potentials fi a — fj, p = /i and a = 1. 

results and the overall picture are in good agreement with 
Refs. 34,35, which were obtained using a different formal- 
ism. 

We have also analyzed the more realistic case in which 
both layers are not equally doped in Fig. 10. Here we 
have tried to emulate the experimental parameters of 
Kim et al. 26 > 27 . We have choosen d = 0.2 and a = 0.2 
which seems to be fitting to their setup. For the dielec- 
tric constant of AI2O3 we have chosen e r w 10. The 
chemical potentials /i a and p p are chosen such that we 
roughly reproduce the variation of densities as shown in 
Fig. 3 of Kim et al. 26 . One can see that the numeri- 
cal result is qualitatively in very good agreement with 
the experiment. In order to make quantitatively accu- 
rate comparison one would have to take into account the 
sample geometry, which we refrain from doing. 

VII. CONCLUSION 

We have made an in depth study of Coulomb drag 
in two parallel monolayers of graphene in a variety of 
regimes ranging from the non-degenerate limit to the 
fully degenerate Fermi liquid limit. On a technical level 
we have employed a description of drag in terms of 



FIG. 10: Numerical modelling of the experiment by Kim et 
al. 2(i with d = 0.2 and a = 0.2. 



the Boltzmann kinetic approach using the variational 
approach within a two-mode approximation, which is 
asymptotically exact in the leading logarithmic approxi- 
mation. We have studied the interplay between interac- 
tions, disorder, and the distance between the two mono- 
layers. Directly at the Dirac point we find the absence 
of drag due to particle-hole symmetry. Still, we find 
that there is an interesting effect of the passive layer on 
the transport properties of the active layer, which comes 
from scattering of electrons and holes in the passive layer 
thereby relaxing a current. This effect is dependent on 
the parameter Td/vF and leads to an interesting temper- 
ature dependence of the inelastic scattering dominated 
single layer conductivity. We point out that this provides 
a promising route towards increasing inelastic scattering 
in graphene bringing closer the collision-dominated hy- 
drodynamic limit. In the close vicinity of the Dirac point 
we found an interesting interplay between elastic and in- 
elastic scattering. We first showed that in the clean limit 
at finite chemical potential there can be a well defined 
finite drag despite the divergence of all individual con- 
ductivities. We find that there is a non-commutativity 
of limits: first taking doping to zero and subsequently 
the disorder yields zero drag pd — while in the re- 
versed order of limits we find p d = —l/a , where cr is 
the single layer conductivity of clean graphene at charge 
neutrality. This effect has been discussed by Schiitt et 
al. 36 in a related work of which we became aware during 
the completion of this manuscript. In the Fermi liquid 
regime, p/T >• 1, we presented a derivation of the drag 
resistivity in a very simplified setting of the Boltzmann 
equation which recovers the standard formula of drag as 
it has been derived in the context of Fermi liquids with 
the Kubo formalism. In the limit T — > we find the stan- 
dard Fermi liquid behavior of the T 2 type, with a distance 
dependence interpolating from l/d° for kpd <C 1 to the 
more l/<i 4 -behavior in the opposite limit. Interestingly, 
we find a previously not discussed regime of T-linear be- 
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havior and 1/d 5 distance dependence if Td/vp > 1. This 
behavior might be relevant for the understanding of re- 
cent experiments on graphene, where deviations from the 
T 2 -behavior have been observed 26 ' 27 . We point out that 
a similar behavior has also been seen in two dimensional 
electron gases and the crossover scales seem compati- 
ble 23 . We closed with a discussion of full crossover curves 
between both regimes which show good qualitative and 
quantitative agreement with experiments if screening is 
taken into account properly. 



as very useful discussions with M. Garst, A. Rosch, M. 
Schiitt, and M. Vojta. We are thankful to A. Geim and L. 
Ponomarenko for making unpublished material available 
to us. This work was supported by the Emmy-Noether 
program FR 2627/3-1 (LF). 



Acknowledgements 



Appendix A: The scattering integral 
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The interaction part of the scattering integral reads 



7g° = 2tt 



where 



d 2 h d 2 k 2 



(2tt) 2 (2tt) 2 1 
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- [1 - / A a (k,t)][l - mk u t)}fZ(k + q.tj/ftk! - q,t)}| 

7 c - 27r 7 (2.) 2 (2.) 2 \ 

<5(i; F fc - v F h - v F \k + q| + v^k, - d\)R n {tf(k, t)^ A (ki, t)[l - tf(k + q, i)][l - ,^ A (k 1 - q, t)] 

[1 - / A a (k,t)][l - /! A (k l5 i)]/ A a (k + q,*)/^ - q,i)} 
+% F fc - « F fci - « F |k + q| + VF |ki - q|)£ 12 {/£(k, i)/! A (k l5 i)[l - /£(k + q, i)] [1 - /» A (k! - q, t)] 

- [1 - / A °(k,t)][l - f_ x (k u t)]f x (k + q.tJZ-^k! - q,t)} 

+% F fc + « F fci - t, F |k + q| - WF |kx - q|)£ 2 {/ A °(k, i)./f (ki, i)[l - / A a (k + q, t)][l - /*(kx - q, t)] 

- [1 - #(k, i)] [1 - f x (k u t)}K(k + q, 0/^ki - q, *)} }, (Al) 



i?i = 4iV|T + __ + (k,k 1 ,q)| 2 and 
+ |T + _+_(k,k 1 ,k 1 -k-q)| 2 
-4T + __ + (k,k 1 ,q)Tj;_ + _(k,k 1 ,k 1 -k-q) 
-4T;_ + _ (k, k x , k x k q)T+__ + (k, ki , q) 

( A2 ) R2 = 4AT|T ++++ (k,k 1 ,q)| 2 

-4T ++++ (k, ki, q)I* +++ (k, k lt ki - k - q) , 
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while 

Ru = 47V|T + __+(k,k 1 ,q)| 2 

R12 = 47V|T + _+_(k,k 1 ,k 1 -k-q)| 2 

R 2 = 47V|f ++++ (k,k 1 ,q)| 2 . 

We have used 
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{K* 2 -Q*)K 2 
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In order to obtain I^F and Iq 1 one simply has to change 
the individual indices. The collision integral due to dis- 
order assumes the form 



(A4) IZ = 2nj ^J(fe-fc 1 )|I7 AA (k,k 



iJI 2 x 



r x (k,t)d m^,t)) - a - /^(k,t))/^(k 1 ,t) 

(A7) 



where again I^f s is obtained from a simple change of in- 
dices. 



Appendix B: The scattering matrix 



We define a space of modes according to 



^(k)=A^Xo( P and ^"(k) = f k 
In this basis the elements of the scattering matrix assume the following forms 
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and 
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C pp and C pa are obtained by a simple exchange of a and p. The full equation can be cast in the form 



D a 





(jaa 

(jpa (jpp 



Xa 



Xa 

X P 



(jaa (jap 
Qpa Qpp 



D a 





(B3) 



(B4) 



1 A. H. Castro Neto, F. Guinea, N. M. R. Peres, 109 (2009). 

K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 



21 



2 X. Du, I. Skachko, F. Ducrr, A. Luican, E. Y. Andrei, 
Nature 462, 192-195 (2009). 

3 K. I. Bolotin, F. Ghahari, M. D. Shulman, H. L. Stormer, 
and P. Kim, Nature 462, 196-199 (2009). 

4 J. Gonzalez, F. Guinea, and M. A. H. Vozmediano, Nucl. 
Phys. B 424, 595 (1994); Phys. Rev. B 59, R2474 (1999). 

5 D. C. Elias, R. V. Gorbachev, A. S. Mayorov, S. V. Mo- 
rozov, A. A. Zhukov, P. Blake, L. A. Ponomarenko, 
I. V. Grigorieva, K. S. Novoselov, F. Guinea, A. K. Geim, 
Nature Phys, 7, 701-704 (2011). 

6 D. E. Sheehy and J. Schmalian, Phys. Rev. Lett. 99, 
226803 (2007). 

7 M. S. Foster and I. L. Aleiner, Phys. Rev. B 77, 195413 
(2008). 

8 A. B. Kashuba, Phys. Rev. B 78, 085415 (2008). 

9 L. Fritz, J. Schmalian, M. Miiller, and S. Sachdev, Phys. 
Rev. B 78, 085416 (2008). 

10 M. Miiller and S. Sachdev, Phys. Rev. B 78, 115419 (2008). 

11 M. Miiller, L. Fritz, and S. Sachdev, Phys. Rev. B 78, 
115406 (2008). 

12 M. Miiller, J. Schmalian, and L. Fritz, Phys. Rev. Lett. 
103, 025301 (2009). 

13 M. Vojta, L. Fritz, R. Bulla, Eur. Phys. Lett. 90, 27006 
(2010). 

14 M.B. Pogrebinskii, Sov. Phys. Semicond. 11, 372 (1977). 

15 P.J. Price, Physica 117B, 750 (1983). 

16 L. Zheng and A.H. MacDonald, Phys. Rev. B 48, 8203 
(1993). 

17 A.-P. Jauho and H. Smith, Phys. Rev. B 47, 4420 (1993). 

18 A. Kamenev and Y. Oreg, Phys. Rev. B 52, 7516 (1995). 

19 K. Flensberg, B.Y.-K. Hu, A.-P. Jauho, and J.M. Kinaret, 
Phys. Rev. B 52, 14761 (1995). 

20 S.M. Badalyan, C.S. Kim, G. Vignale, and G. Senatore, 
Phys. Rev. B 75, 125321 (2007). 

21 R. Asgari, B. Tanatar, and B. Davoudi, Phys. Rev. B 77, 
115301 (2008). 

22 For a review see A.G. Rojo, J. Phys.: Condens. Matter 11, 
R31 (1999). 

23 T.J. Gramila, J. P. Eisenstein, A.H. MacDonald, L.N. Pfeif- 
fer, and K.W. West, Phys. Rev. Lett. 66, 1216 (1991). 

24 U. Sivan, P.M. Solomon, and H. Shtrikman, Phys. Rev. 
Lett. 68, 1196 (1992). 



25 J. Moore, Nature Phys. 5, 378 (2009); J.E. Moore, Nature 
464, 194 (2010); M.Z. Hasan and C.L. Kane, Rev. Mod. 



Phys. 82, 3045 (2010); X.-L. Qi and S.-C. Zhang, ibid. 83, 
1057 (2011). 

26 S. Kim, I. Jo, J. Nah, Z. Yao, S.K. Banerjee, and E. Tutuc, 
Phys. Rev. B 83, 161401 (R) (2011). 

27 S. Kim and E. Tutuc, Solid State Comm. (2012), 
http://dx.doi.Org/10.1016/j.ssc.2012.04.032. 



R. V. Gorbatchev et ai, to appear. 



28 

29 W.-K. Tse, Ben Y.-K. Hu, and S. Das Sarma, Phys. Rev. 
B 76, 081401(R) (2007). 

30 B.N. Narozhny, Phys. Rev. B 76, 153409 (2007). 

31 M.I. Katsnelson, Phys. Rev. B 84, 041407(R) (2011). 

32 N.M.R. Peres, J. M.B. Lopes dos Santos, and A.H. Castro 
Neto, Europhys. Lett. 95 18001 (2011). 

33 E.H. Hwang, R. Sensarma, and S. Das Sarma, Phys. Rev. 
B 84, 245441 (2011). 

34 B.N. Narozhny, M. Titov, I.V. Gornyi, and P.M. Ostro- 
vsky, arXiv: 11 10.6359. 

35 M. Carrega, T. Tudorovskiy, A. Principi, M. I. Katsnelson, 
and M. Polini, arXiv:1203.3386 (2012). 

36 M. Schiitt, P. M. Ostrovsky, M. Titov, I. V. Gornyi, 
B. N. Narozhny, and A. D. Mirlin, arXiv:1205.5018 (2012). 

37 J. C. W. Song and L. S. Levitov, arXiv:1205.5257 (2012). 

38 C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sor- 
genfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shep- 
ard, and J. Hone, Nature Nanotech. 5, 722-726 (2010). 

39 J. M. Ziman, Electrons and Phonons, Oxford University 
Press, Oxford (1960), Chapter 7. 

40 M. Schuett, P. M. Ostrovsky, I. V. Gornyi, and A. D. Mir- 
lin, Phys. Rev. B 83, 155441 (2011). 

41 B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New J. 
Phys. 8, 318 (2006). 

42 E.H. Hwang and S. Das Sarma, Phys. Rev. B 75, 205418 
(2007). 

43 L. Fritz, Phys. Rev. B 83, 035125 (2011). 

44 A. Larkin and A. Varlamov, Theory of Fluctuations in Su- 
perconductors (Clarendon Press, Oxford, 2004). 

45 L. Ponomarenko, private communication. 

46 A. Levchenko and A. Kamenev, Phys. Rev. Lett. 101, 
216806 (2008). 



